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ABSTRACT 


This paper describes the basic navigation solution for position and 
velocity based on range and delta range (Doppler) measurements from NAVSTAR 
Global Positioning System satellites. 

The application of discrete filtering techniques is examined to 
reduce the white noise distortions on the sequential range measurements. A 
second order (position and velocity states) Kalman filter is implemented to 
obtain smoothed estimates of range by filtering the dynamics of the signal 
from each satellite separately. 

Test results using a simulated GPS receiver show a steady-state noise 
reduction, the input noise variance divided by the output noise variance, 
of a factor of four. 

Recommendations for further noise reduction based on higher order 
Kalman filters or additional delta range measurements are included. 
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INTRODUCTION TO GPS NAVSTAR 


I . 


The Global Positioning System (GPS) NAVSTAR is a satellite based 
navigation system. The system provides the user with highly precise 3D 
position and velocity information and time worldwide on a continuous, all- 
weather basis. The system consists of a Space Segment, a Control Segment 
and a User Segment. 

When fully operational in 1990, the Space Segment will consist of 18 
satellites in six 12 hour orbits of 3 satellites. Each satellite will con- 
tinuously broadcast a message containing precise information relative to 
its own position (ephemeris) and clock accuracy and less precise infor- 
mation relative to the entire constellation position (almanac). 

The Control Segment consists of monitor stations and a master control 
station. The monitor stations transmit satellite tracking data to the 
master control station, which determines the satellites' orbital parameters 
and communicates them to the satellites for retransmission to the users. 

The User Segment consists of the equipment necessary to derive posi- 
tion, velocity and time from the information received from the satellites. 

Four satellites are normally required for navigation purposes, and 
the four offering the best geometry can be selected manually or automati- 
cally by receivers using ephemeris information transmitted by the satel- 
lites. Ranges to the four satellites are determined by scaling the signal 
transmit time by the speed of light. The transmitted message contains 
ephemeris parameters that enable the user to calculate the position of each 
satellite at the time of transmission of the signal. 

Operation of the system requires precise synchronization of space 
vehicle (SV) clocks with GPS time, which is accomplished by the use of an 
atomic frequency standard in each space vehicle and the use of correction 
parameters that are provided by the Control Segment. The requirement for 
users to be equipped with precision clocks is eliminated by the use of 
range measurements from four satellites. If users maintained precision 
clocks synchronized with GPS time, navigation could be accomplished with 
only three satellites. In that case, the user could be thought of as being 
at the intersection of three spheres, with centers located at the satel- 
lites. The fourth satellite permits an estimate of the user's clock error. 
In this case, the user position contains four unknowns consisting of posi- 
tion in three dimensions and the error, or fixed bias, in the user's impre- 
cise clock, which can be solved by simultaneous solution of the four 
equations . 

The range measurements from satellite to user are corrupted by 
several distortions. Basically, these distortions can be divided into two 
groups; random and bias distortions. The random distortions can be reduced 
by use of filtering techniques (see Chapters 4, 6 and 7). It is not 
possible to make a distinction between the actual range and the bias 
errors. Therefore, bias distortions can not be reduced by use of filtering 
techniques. 


- 1 - 



However, there is a way to obtain a very accurate estimate of the 
bias distortions, by implementing the so-called differential GPS (DGPS) 
approach . 

The DGPS approach uses a receiver/transmitter at a known, fixed loca- 
tion. This receiver compares its GPS derived position to its actual sur- 
veyed location, and transmits the errors to suitably equipped users to 
allow them to improve their own solutions. Various DGPS concepts can be 
implemented [1], 

In this paper the accent is on the following concept: 

- a receiver is placed in a known location and the errors R(err) 

(bias errors) from user to satellite are measured: 

R(err) = (measured range) - (true range) 

The correction term, R(err), is then transmitted to the user 
receiver . 
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II. GPS SIGNAL STRUCTURE 


The GPS signal consists of two components, Link 1 (LI) at a center 
frequency of 1575.42 MHz and Link 2 (L2) at a center frequency of 1227.6 
MHz. Each of these two signals LI and L2 , is modulated by either or both a 
10.23 MHz clock rate precision P signal and/or by a 1.023 MHz coarse/ 

acquisition C/A signal. Each of these two binary signals is formed by a 

P-code or a C/A code which is modulo-2 added to 50 bps data D, to form P e 

D and C/A © D. The PSD and C/A ® D signals are modulo-2 added to LI and 

L2 in phase quadrature. 

The P-code is a pseudo-random sequence with a period of 1 week. The 
C/A code is a unique Gold code with a period of 1 msec. The user has the 
capability to duplicate both P and C/A codes and the transmission time is 
determined by measuring the offset that has to be applied to the locally 
generated code to synchronize it with the code received from the satellite. 

Since the P-code has a wider bandwidth, it will be more difficult to 
acquire than the C/A code, but it will provide greater accuracy and addi- 
tional anti-jam protection. 

In this report, we are only concerned with the C/A code, as the P- 
code will only be available for authorized users. However, the C/A code, 
in combination with differential GPS, will provide an accuracy on the same 
order of magnitude as the P-code. 
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III. DEFINITIONS (PSEUDO) RANGE AND DELTA (PSEUDO) RANGE 


The principle of GPS navigation is based on the delay between the 
time the navigation code is transmitted by the satellite and the time the 
code is received by the user. This time delay can be considered as a range 
from satellite to user in a vacuum, denoted as the true range: 


R = c * ( Tr - Tt ) meters 

where c « GPS value for speed of light (2.99792458 * 10**8 m/s) 

Tr= GPS receive time 
Tt* GPS transmit time 

However, there are additional time delays, including: 

- the user clock bias TBu with respect to true GPS time; 

- atmospheric delay TDa, as the signal is not transmitted in a 
vacuum ; 

- receiver delay TDr between the antenna phase center and the 
code-correlation point in the receiver; 

- the transmit clock bias TBs with respect to true GPS time 

The range measured by the user is called pseudorange and is defined as: 

PR = R + c * ( TDa + TDr + TBu - TBs ) meters 

The change in true range over a specified time interval is called delta 
range : 

DR = R(t + dt) - R(t) 

where dt* the time interval between the two range measurements 

Delta range divided by the specified time interval is the mean rela- 
tive velocity between the user and satellite for that time interval. 

The delta pseudorange is defined as the change in pseudorange over a 
specified time interval: 

DPR = PR ( t + dt) - PR(t ) 

= DR - c * ( DTBs - DTBu) 

where DTBs * the satellite clock bias change over the time interval dt 
DTBu * the user clock bias change over the time interval dt 

Delta range is usually obtained by counting the number of carrier 
cycles that occur over a finite interval in the receiver phase lock loop. 
The total phase change (Doppler count) accumulated over the measurement 
interval is equivalent to the integral of the Doppler phase rate 
(frequency) due to relative motion between the user and the satellite 
during the measurement interval . 
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t 2 


Dopplercount AN = 


f g - f^ 1 ) 


dt 


(3.1) 


tl 


where t£ = tj + At 

At = Dopplercount interval 


fg = precise ground station 
frequency 

f r (t) = received frequency 


The mean velocity during the measurement interval is given by: 
AN*A 


v = 


At 


(3.2) 


where A = carrier wavelength 
and the delta pseudorange is given by: 
DPR = AN * A 


(3.3) 
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IV. GPS ERRORS 


The GPS navigation accuracy is a function of the product of two com- 
ponents : 

- Measurement Errors 

- Geometric Dilution of Precision GDOP 

A. Measurement Errors 

Error contributions have been attributed to the various system 
segment contributors; the Space segment, the Propagation Link and the User 
segment. The user equivalent range error (UERE ) and the user equivalent 
delta range error (UEDRE) are measures to describe these contributions to 
pseudorange and delta pseudorange measurements. 

B. Geometric Dilution of Precision GDOP 

GDOP is the degradation of accuracy that results when the geometry of 
the satellites is not optimal: 

GDOP = J o x 2 + q y 2 + (7 Z 2 + cr t 2 / 0p 


where a x 2 , <jy z , a 2 2 are the variances of user position 
^t 2 is the variance of user time 
Op is the pseudorange standard deviation 

GDOP can be partitioned into separate position and time variances: 

GDOP 2 = PDOP 2 + TDOP 2 

where PDOP = J a x 2 + a y 2 + o z 2 / a p 

= Position Dilution of Precision 
TDOP = a t /<7p 

= Time Dilution of Precision 

In a navigation system PDOP is used for the determination of the 
position accuracy, knowing time is generally not necessary. PDOP is deter- 
mined by the satellite constellation, user's geographical location, mask 
angle, and time of the day [2]. A typical value of PDOP is 3.4 with a 
probability of 0.9 where the elevation mask is 5 degrees for a 24-satellite 
constellation . 
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C. Position Error 


The rms position error (x,y,z) 



is related to the rms radial range error <7 R by 
a p = PDOP * a R 

where a p = 1 a position error in three dimensions 

a R = square root of sum of squares of UERE contributions 
PDOP = position dilution of precision in three dimensions 

D. Velocity Error 
The rms velocity error (x,y,z) 



is related to the rms radial delta range error ojjp by 

a v = PDOP * o‘dr 

where a v = 1 a velocity error in three dimensions 

a DR = square root of sum of squares of UEDRE contributions 

E. Influence of UERE and UEDRE on Range and Delta Range 

The influence of the UERE and UEDRE on the range and delta range 
(Doppler) can be divided accordingly to the summary given in table 1 [3-6]. 

Using a differential GPS approach, it follows from table 1 that the 
expected errors in range (position) and delta range (velocity) due to ran- 
dom noise are respectively 14 meters and 1.3 centimeters. 
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Table 1. Influence of the UERE and UEDRE on the Range and 
Delta Range (Doppler) for the C/A Code. 


ERROR CONTRIBUTORS 

RMS RANGE 
ERROR 
IN M 

RMS DELTA 
RANGE 

ERROR IN M 

STATISTICS 

COMMENTS 

satellite ephemeris 

1.5 


bias 

uncorrelated 
between SV's 

satellite group and 
clock relativistic 
frequency shifts 

1.0 


bias 

uncorrelated 
between SV's 

satellite clock noise 


0.003 

white noise 

cesium clock 

receiver noise and 
resolution 


0.008 

white noise 

reference [5] 

receiver noise 

< 10 ** 


markov 

reference [5] 
and [4] 

receiver resolution 

0.0005 


white noise 

reference [5] 

ionospheric and 
tropospheric delays 

< 25 


bias 


propagation gradient 


0.01 

markov 


multipath error 

10 


white noise 



** This error consists mainly of the thermal tracking jitter error, which 
decreases when C/No increases. In case an accurate measure of C/No in the 
receiver and the bandwidths is available, a prediction of this error in 
real time can be obtained [5] . 
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V. THE NAVIGATION SOLUTION 

A. Introduction 

This Chapter gives an overview of the basic Methods to derive user 
position, velocity, clock bias and clock bias rate from pseudorange and 
delta pseudorange (Doppler) measurements. Some of these algorithms were 
implemented in the computer language FORTRAN IV and tested on a IBM 
System/370 in combination with a GPS simulation program. 

In order to simulate the GPS using an 18-satellite constellation, two 
computer programs were used [7]: 

- the NAVSTR program which provides satellite position, satellite 
velocity, pseudorange and delta pseudorange. The program has the 
option to introduce several distortions into the pseudorange data. 

- the USRPOS program which simulates a helicopter flight pattern con- 
taining several maneuvers and provides user position and velocity. 

The simulations used in this paper are all based on a flight duration 
of 270 seconds. The time interval between fixes is 0.3 second. During 
the first 138 seconds the USRPOS program simulates an unaccelerated 
straight level flight. During the next 31 seconds a turn and a -2 
knots/hour acceleration is simulated. The last part is again a straight 
level flight. 

All the programs are written in double precision (precision of 
approximately 16.8 digits [8]), instead of single precision (precision of 
7.2 digits) because of the large range and delta range values with respect 
to the desired computation precision. For instance, a delta range value of 
9000 nmi requires at the very least a computation accuracy of 8 decimals to 
obtain 1 meter computation precision. 

B. Position 

The GPS navigation position fix and time bias can be obtained from 
the following basic equations (in Earth Centered Earth Fixed (ECEF) 
coordinates ) : 


(r 2 - b) 2 = 
(r 2 - b) 2 = 
(r 3 - b) 2 = 


< u x - s xi ) 2 
< u x " s x2 ) 2 
< u x " s X 3) 2 
( u x ” s x4 ) 2 


+ (°y - Syl) 2 

+ (Uy - Sy 2 ) 2 

+ (Uy - Sy 3 ) 2 

+ (Uy ~ S y 3 ) 2 


+ (u z - s zl ) 2 
+ (u z - s z2 ) 2 

+ (u z - S z3 ) 2 

+ (u z - s z4 ) 2 


(5.1) 


(r 4 - b) 2 



where (u x , u y , u z ) is the user position; b is the user clock bias with 
respect to GPS tine; rj is the true range distance between the user and ith 
satellite; (s x j, s y j , s 2 |) is the position of the ith satellite. 

The user state (u x , u y , u z , b) can be determined by solving eq. 5.1 
directly, given the measured pseudoranges R ■ (rj r 2 r 3 r 4 ) T and the three 
satellite position components (s x j, s y i , s z i), i = 1, 2, 3, 4. This 
approach is computationally unwieldy, especially because the calculations 
are usually executed in a coordinate system not very close to the problem. 

A good alternative approach is to linearize (5.1) about an estimate of the 
user state and solve successively for position corrections based on new 
measurements . 

The basic navigation equations can be linearized by employing incre- 
mental relationships as described by Jorgensen [9] . The resulting 
linearized navigation equations are given by: 


Un7 ~ S' 


u nx s xi . u ny s yi u nz “ a zi . .. . 

Ax + Ay + Az + Ab = Ar, 


m 


- b. 


ni 


- b, 


r ni bn 


(5.2) 


where u nx , u ny , u nz , b n are nominal (a priori best estimate) values of 
u x , u y , u z , b; Ax, Ay, Az, Ab are the corrections to these nominal values; 
r ni are the nominal pseudorange measurements from the ith satellite 
(i = l ,2,3,4) ; Ar^ is the difference between the actual and nominal measure- 
ments . 

The linearized equations can be expressed in matrix notation as follows: 


H-AU = AR or AU = IT 1 AR 
where AU = (Ax Ay Az Ab) T 

AR = (Arj Ar 2 Ar 3 Ar 4 ) T 
H = (h a h 2 h 3 h 4 ) T 

[ u nx " s xi u ny ~ s yi u nz “ s zi 

h i = * 

L r ni - b n b n r n i - b n 


(5.3) 


(5.4) 


The known quantities Arj are actually incremental pseudorange 
measurements. They are the differences between the actual measured 
pseudoranges and the measurements that had been predicted by the user's 
computer based on the knowledge of satellite position and the user's most 
current state. The quantities to be computed, Ax, Ay, Az and Ab are 
corrections that the user will make to his current state. 
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C. Position Algorithm 


The methodology for calculating the user state based on the linearized 
equations (5.3) proceeds as follows [10]: 

1) Guess initial user state U n = (u nx , u ny , u nz , b n ) 

2) Obtain rj and s x j , s y i , s Z j for i = 1,2 ,3, 4 

3) Compute r n < = /, ~ , 

/ ( u nx -s xi) 2 + ( u ny ~ s yi) 2 + ( u nz ~ s zi) 2 + b n 

4) Compute Ar* = rj - r n j for i = 1,2, 3, 4 

5) Compute H from (5.4) 

6) Compute G = H -1 

7) Compute AU = G * AR 


8) Update state estimate u n = ^n + AU 

9) if | AU | > e, a specified tolerance, go to step 3) for another 
iteration, else go to step 2) for the next position + clock bias 
fix 

The program has a build-in limit for the maximum number of iterations 
per fix. 

Figure 1 shows the results for the 3-dimensional position error where 
e = 0.001, and no pseudorange distortions were introduced. The number of 
iterations per fix to reach the desired accuracy during this simulation was 
one or two. 

Figure 2 shows the magnitude of the 3-dimensional position error when 
multipath error was introduced. White noise statistics are assumed for the 
multipath error having a zero mean and a variance of 10 meters. 

D. Velocity 

The user velocity can be derived from the range information 
(position) or from the carrier Doppler frequency. The Doppler velocity is 
expected to have a much greater accuracy (approximately 1.3 centimeters) 
than the velocity derived from range information. However, to obtain the 
Doppler velocity, additional hardware is required. 
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1. Velocity Derived from Carrier Doppler Frequency 

Based on the velocity equation 3.2 we can derive a similar algorithm 
as was given for the position fix. The user state consists of the user 
velocity (v x , Vy, v z ) and the clock bias rate b. The velocity measured on 
the 1 ine-of-sight from satellite to user is called r. In the simulation, 
the satellite velocity is provided by the NAVSTR program. In an actual 
receiver, the satellite velocity can be easily derived from satellite posi- 
tion within the desired accuracy [ 11 ]. 


The algorithm for the user velocity computation, based on a similar 
linearization about the user state as used for the position fix, proceeds 
as follows: 


1) Guess initial user state U n = (v nx , v ny , v nz , b n ) 

2) 

3) 


+ b n 


where rangeflag is the sign of the change in range 


4) Compute Ar^ = ri - r n j for i = 1,2, 3, 4 


5) Compute H = (hj h 2 h 3 I 14 ) T 

v nx " sv xi v ny - sv yi v nz " sv zi 
where h^ » 1 

r i " b n ri - b n r i “ ^n 


6 ) Compute G = H _1 


7) Compute AU = G * AR 


where AR = (Arj Ar 2 Ar 3 Ar 4 ) T 
AU = (Av x AVy Av z Ab) T 


Obtain r^ , sv xi , sv yi , sv zi for i = 1,2, 3 , 4 

Compute r ni = (rangeflag * f. ~~ ' ) + 

/ ( sv xi“ v nx ) 2 + (sv yi -v ny )2 + (sv zi -v nz )2 } 
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8) Update state estimate U n = U n + AU 

9) If | AU | > c, a specified tolerance, go to step 3) for another 
iteration, else go to step 2) for the next velocity + clock bias 
rate fix. 


The performance of the velocity Doppler algorithm is similar to the perfor- 
mance of the position algorithm as expected (see figure 1). 

2. Velocity derived from range information 

The basic method to derive the velocity from range information is to 
fit a nth order function through n+1 position measurements based on range 
information. 

2.1. First-order Estimation 

The estimated velocity over the time interval between two successive 
position fixes equals: 


x(T + nT) - x(nT) 

v - 

T 


The first- order estimation provides a reasonable estimate of the user 
velocity, however the velocity estimates are not available at the same time 
as the position fixes. 

2.2. Second-order Estimation 

In fact, this method supposes a constant acceleration for the 
(vehicle + satellite) dynamics by fitting a second-order function 
(parabola) through three successive position fixes. The parabola provides 
an estimation of the position at some intermediate time, t, within the 
measurements : 


y(t) * aj + a 2 t + a 3 t 2 te [nT, nT + 2T) 

Given the three measurements: 

yi (nT) = aj + a 2 nT + a 3 (nT) 2 

y 2 (nT + T) = a 1 + a 2 (nT + T) + a 3 (nT + T) 2 

y 3 (nT + 2T ) = a 1 + a 2 (nT + 2T) + a 3 (nT + 2T) 2 
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it can be easily be shown that 


a l 


yi (2 + 3n + n 2 ) + y 2 (-4n - 2n 2 ) + y 3 (n + n 2 ) 

2 


a 2 


y 1 {-3 - 2n) + y 2 (4 + 4n) + y 3 (-2n -1) 

2T 


yi - 2 y 2 + y3 

a 3 » 

2T 2 


The velocity estimate at some intermediate time te [nT, nT + 2T] is given 
by 


v ( t ) = a 2 t + 2a 3 t 


Using a second-order approximation, acceptable velocity estimations 
are obtained as long as the range data are not corrupted by distortions. 
When introducing distortions, the velocity errors will be in the same order 
of magnitude or larger than the position errors for a time interval between 
fixes smaller than 1 second. 
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VI. GPS SIGNAL FILTERING 


ORIGINAL 


OF POOR 


QUALITY 


DGPS is capable of removing the bias component in the range data. 

The random distortion however will still degrade the range data. For 
example see figure 3 for a multipath distortion of the range data. 

There exist two extreme methods [12] for filtering the signals from 
the satellites and forming the navigation solution for position, velocity, 
perhaps acceleration, and clock parameters. The most powerful method is to 
input the raw observables of pseudorange and delta pseudorange from all 
satellites observed into a total solution filter which solves simulta- 
neously for all the navigation parameters. The filter is typically a 
Kalman filter, or some modification thereof. 

The second extreme method is to filter the dynamics of the signal 
from each satellite separately, forming smoothed estimates of pseudorange, 
velocity, and perhaps acceleration, and then form a memoryless flash solu- 
tion for the navigation parameters, without further smoothing. This method 
is not as powerful as a general Kalman filter, but often has distinct 
implementation advantages without significant loss in performance. Many 
navigation receivers use a combined method, in which the individual signals 
are tracked and smoothed in phase locked loops , and the loop outputs are 
further smoothed in the navigation filter. This combined method has the 
advantage over the first method that the filter update rate does not have 
to be as high because the raw data are smoothed. 

In this paper we will only consider discrete filtering techniques, as 
the range data are a discrete measurement sequence. The best way to imple- 
ment discrete filters is by using a digital computer. Some important con- 
siderations in using this approach are: 

- computer time: for instance the time required for matrix operations 
in a Kalman filter grows as the cube of the state dimension [13]; 
also a higher calculation precision requires more computer time. 

- computer storage: for instance the storage required for the Kalman 
covariance matrix grows as the square of the state dimension [13]. 

For the purpose of establishing feasibility and a performance bound, 
it is sufficient to consider the second filter type: separate tracking 
filters for each satellite, followed by a flash navigation solution as pre- 
sented in Chapter 5. 

It is also sufficient to consider the tracking filter only, since the 
maximum navigation filter position error due to dynamics is equal to the 
maximum position error in one tracking filter, and since the error in posi- 
tion solution due to random noise is equal to the pseudorange error 
multiplied by the PDOP which does not depend on filter type [12]. 

A. One Channel Filtering 

Various techniques are available for the separation of the range 
signal from the random distortions [14] and [15] . In this paper we will 
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mainly consider a particular branch of filters, the so-called least-squares 
filters. Simply stated, the least-squares filter problem is this: given 
the spectral characteristics of an additive combination of signal and 
noise, what linear operation on this input combination will yield the best 
separation from the noise, where best means minimum mean-square error [14]. 
Besides the least-squares filtering one other type of filter will be men- 
tioned, the second-order filter as presented by Hurd [12]. 

1. Second-order Filter 

W. J. Hurd [12] describes a second order filter with both poles at 
the same location, using only range and not delta range as the input. This 
filter is characterized by one parameter, and when the update rate is suf- 
ficiently high this parameter is effectively a time constant T. This is 
valid when the time constant between updates is small compared to T. The 
response of the filter to random noise is to reduce the rms error by a fac- 
tor of 1.12 times the square root of (T divided by the time between 
inputs ) . 

In case of an update time of 0.3 second and a time constant of 2 
seconds, the raw range errors are reduced by a factor of 2.6. However, the 
position error due to acceleration is equal to acceleration, (a), times T 
squared. A time constant of 2 seconds will result in a range error due to 
acceleration of 4a meters. 

2. Wiener Filter 

Reference [14] gives a very good description of the Wiener filtering 
technique. The main disadvantages using a discrete Wiener filter are: 

- the filter is non-recursive, so there is a 'growing memory' problem 
caused by the need to store all of the past measurement data. This 
is very unwieldy in on-line applications. 

- the filter is difficult to extend to more complicated multiple- 
input multiple-output problems. 

3. Kalman Filter 

The main features of the discrete Kalman filter formulation and solu- 
tion of the filtering problem are [14]: 

- vector modeling of the random processes under consideration; 

- recursive processing of the noisy measurement (input) data 

Beside these features, the Kalman filter approach allows many exten- 
sions and variations for specific implementations. 

B. One Channel Kalman Filtering 

The discrete Kalman filter is for the range filtering problem in this 
paper superior, to both the Wiener and the second order filter. The Wiener 



filter is having the 'growing memory' problem and both the second-order 
filter and the Wiener filter are difficult to extend to a filter having 
range and delta range as input data. 

The approach in which the very accurate delta range derived from 
Doppler measurements is also fed into the filter is expected to reduce the 
remaining random noise error in the position compared to the filter with 
range measurements only. 
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VII. DEVELOPMENT AND IMPLEMENTATION OF A DISCRETE KALMAN FILTER 
A. Description of the Discrete Kalman Filter 

The derivation of the discrete Kalman filter is very well documented 
[14 - 16], therefore in this paper we will only present the final filter 
equations and the necessary definitions. 

If the true state of a Markov process and a measurement are defined . 
by the equations: 


£ k + 

1 = h?k + Sk 

(7.1) 

2 k = 

H k x r + v k 

(7.2) 


where Xk is the state vector at time tj* 

is the state transition matrix, representing the known 
system dynamics 

Zk is the measurement vector (filter input) at time tj< 

Hr is the measurement matrix, representing the relationship 
between the measurement zr and the state vector xr in 
the absence of noise 

wr is the plant noise with covariance matrix Q 
Yk is the measurement noise with covariance matrix R 
E l *k ¥i T ] = 0 for all k and i 

then the best estimate xr the system state vector and its covariance 
matrix Pr can be updated by the equations: 

x k+1 (-) * <f> k + f 

P k +l*-> ” * k p k (+ * * k T + Q k (7.3) 

P k < + ) = (I - K k H k ) P k <-> 
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Figure 4. Coaputational Steps. 
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Figure 5. Block Diagram of System Implementation. 




^ + Kk (zk - Hk x ^ ) 


• 4) 


K k = p k M H k T { H k Pk<-)H k T + R k I" 1 


where Xk^ > p k^ and 2Sk^ + ^» p k^ + ^ are respectively estimates of Xk and 
Pk just before and just after tine tk- 
Kk is the Kalman gain (blending factor) 

Figure 4 shows the shows the sequence of the computational steps pic- 
torially and figure 5 shows a block diagram of the system model, the 
measurement and the discrete Kalman filter. 

B. Kalman Filter Implementation Considerations 

When the discrete Kalman filter equations (7.1 - 7.3) are applied to 
practical problems, several difficulties quickly become obvious: 

- the optimal filter must model all error sources in the system 
under consideration; 

- it is assumed in the filter equations that exact descriptions 
of the system dynamics, error statistics, and the measurement 
process are known. 

As an unlimited computer capability is not usually available most of 
the implemented Kalman filters are sub-optimal; less than the optimal 
number of states are included in the system and/or measurement model and no 
exact description of the error statistics is available. The performance of 
a sub-optimal Kalman filter is not yet based on a unified body of theory 
and practice. Fortunately there are many examples of discrete Kalman 
filter implementations available (14 - 16] from which concepts for a suc- 
cessful implementation can be derived. 

Some of these concepts are: 

- derive a system model representing the real system as close 
as possible; 

- use at least double precision in the computation of the filter 
equations (especially in off-line implementations) to prevent 

the filter to diverge due to calculation errors (for a description 
of the performance degradation in digitally implemented Kalman 
filters [17]); 

- implement 'reasonableness checks' in the filter, e.g. remove 
unlikely input data, when a matrix is expected to be symmetric, 
force it to be symmetric during or after each operation on that 
matrix ; 

- find an acceptable evaluation of the plant covariance matrix Q. 

In most implementations it is not a trivial task to derive the Q 
matrix from its definition: 


OF POOR 
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Qk = E [ wk Wk T ] 


(7.4) 


tk+l tk+l 


= E { 


(^k+l »u)G(u)E[w(u)w t (v)] * 


*k 


*k 


*GT( v )<|) T (tk+i ,v) dudv} 


where E [w(u)w T (v)] is a matrix of Dirac delta functions that, 
presumably, is known from the continuous model. 

G represents the influence of the noise vector w on the 
measurements . 

C. System Modeling 

In the area of pulse Doppler surveillance radars several discrete 
Kalman filters were implemented [18]. The signal model of the pulse 
Doppler radar is very similar to the GPS signal model, therefore, most of 
the techniques derived for the pulse Doppler radar can be used for the GPS 
one-channel filtering. 

Some of the system models as presented by Fitzgerald [18] are: 

- three state exponentially correlated acceleration (ECA) model; 

- two state exponentially correlated velocity (ECV) model; 

- two or three state integrated white noise or 'random walk' 
acceleration (RWA) model. 

In order to obtain a performance bound and a basis for a feasibility 
study, a two state RWA model is implemented. The RWA tracking problem is 
described by the state-space model: 


*k+l = ‘tec + "k (7.5) 

where xj< = ( r r ) T is the state-vector at time t^ 
r is the range from satellite to user 

r is the relative velocity on the line-of-sight between 

satellite and user 

is a stationary white noise process representing the random 
acceleration noise at time t^ 
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<f>k is the state-transition matrix at time t k 

D. Measurement Modeling 

Two types of filter input data are available: range and delta range. 
In a DGPS approach, the range data are degraded by random noise with an rms 
value of approximately 14 meters (see Chapter 4) and the delta range 
derived from Doppler is accurate up to 1.3 centimeters. In order to com- 
pare the accuracy of the velocity derived from range with Doppler velocity 
and because of defining a basis for a feasibility study, the filter is fed 
with range data only. However, the filter design makes allowance for an 
extension to both range and delta range measurements. 

The measurement equation based on range data only becomes: 


z k = »k5k + v k 


where z k (a scalar) is the range measurement at time t k 
is the measurement matrix at time t k 
v k scalar) is a stationary white noise process representing 

the random measurement noise at time t k 

The measurement matrix H k is independent of the time. In the absence 
of measurement noise, the measurement z k equals the range r, therefore the 
H matrix becomes (10). 

E. Computation of the State Transition Matrix <f> 

The dynamical system as described by equation (7.1) has fixed parame- 
ters; i.e., the continuous state transition matrix P is a constant. In 
this case, the state transition matrix <p may be written as an exponential 
series : 

4>k = e FT = I + FT + + 


where T = t k+1 - t k 
In the continuous case: 


r 1 
1 = 

0 

1 


r 

. r 1 

0 

0 


. r . 


therefore 
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' 1 0 ] [ o t 1 fool [it 

0 = + 4 - +••• = 

0 1 J [ 0 0 J l 0 0 j [ 0 1 , 

F. Computation of the Error Covariance Matrix Q 
The Q matrix is defined by equation (7.4): 

Q = E [wkW}< T ] 

where w^ = [wj W£] T 
T 

Wj = u*a(u)du 

0 

T 

Wg = a(u)du 

0 

a(u) is the random white noise acceleration 
T T 

Qll - u*v* E |a(u)a(v) }dudv 

0 0 

E{a(u)a(v)| = a a 2 T 6(u - v) 

where a a 2 is the variance of the acceleration and 6(u-v) is the Dirac 
delta function 

6 ( u-v ) = 1 if u = v 

6(u-v) = 0 if u * v 

it follows that 
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also 


and 


also 


ment 


Qll = 


T 

1 

v 2 a a 2 T dv = 3 o a 2 T 4 

0 


Q 12 = 


T T 


0 0 


u»E (a(u)a(v) }dudv 


va a «.Tdv = 2 a a 2 T3 


T 

E{a(u)a(v) (dudv 

0 

• T dv = a a 2 T 2 

G. The Measurement Covariance Matrix R 

As the measurement vector zj^ is a scalar, the covariance matrix R is 
a scalar, given by the variance <7 X 2 of the measurement noise. 

H. Expected Filter Performance 

1. Sampling Rate Versus rms Acceleration 

In order to keep the accuracy between observations below the measure- 
(sensor) error, the following equation should be satisfied, from [19]: 


Q22 = 


<V 
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CT a T 2 < 0.24 a x 


where <7 X is the measurement accuracy and 
<J a is the rms acceleration 

The maximum error will occur just before a measurement. In the 
model, the time-constant T is set to 0.3 second, therefore 

/ °x >0.6 
* °a 


For instance, if the measurement accuracy a x = 1 foot, then the accuracy 
between observations will be less than 1 foot if the rms acceleration is 
below 2.78 ft/s*s. 

2. Velocity Accuracy 

However, no direct velocity measurements are made; velocity can be 
derived from position by subtracting two successive position measurements 
divided by the time interval between them. 

The corresponding accuracy in velocity is given by [19]: 
a v ~ ^ p v - J 1 (/“ITS? + d^T 2 


4a x 

where r = 

a a T 2 


For example, if a x = 1 foot; a a = 2.78 ft/s*s; T = 0.3 s, then it follows 
from the above equation that the rms velocity error < 7 V = 1.57 ft/s. 

I. Filter Summary 


The discrete Kalman filter for one channel range data tracking pro- 
ceeds as given by the diagram in figure 4, where: 


*k <_) 


estimated initial range 
estimated initial velocity 
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p k ( -) = 


. P2 (_) P3 (_) 


Pjf - ). P 2 ^^, P 3 ^~^ are initial guesses having 
approximately the same order of magnitude 
as the measurement noise variance. 


H = (10) 


T 4 /3 

TV2 


1 

T 

T 3 /2 

T* 

<P = 

0 

1 


R = a x * 


J. Simulation Results 

The measurement noise is generated by the NAVSTAR simulation program. 
The simulation was used in the mode where multipath error is added to the 
range data. 

Figure 3 shows the raw range data produced by the NAVSTAR simulation. 
The variance of the white noise equals 10 meter [7]. 

Figures 6, 7 and 8 show the filter response for different values of 
the noise variances a x 2 and a a 2 , as summarized in table 2. 

The calculated parameters for optimal filter response are given by 
the first row of table 2. 


Table 2. Single Channel Discrete Kalman Filter Performance 
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Figure 6. Filter Response; o x z 
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Figure 8. Filter Response; 



VIII. CONCLUSION AND RECOMMENDATIONS FOR FURTHER RESEARCH 

The discrete Kalman filter for tracking GPS signals separately is a 
very flexible method for an adequate solving of the random noise problem. 

Although a very simple filter (range measurements only) was used in 
this paper, the (steady state) noise variance was reduced with a factor 4 
(table 2 ) . 

In case the highly accurate delta range (Doppler) measurements are 
available to the tracking filter, the resulting improvement in position 
estimation can be an order of magnitude or more [20]. The filter is also 
relatively easy to extend to a three state filter capable of handling the 
Markov noise produced by the receiver (see table 1) as well. Another 
possibility is to add additional states for allowing more complex vehicle 
dynamics [18, 21, 22]. 

The results shown in figures 6, 7 and 8 are also representative for 
the final three dimensional position error; see Chapter 6, also compare 
figure 2 (magnitude of the 3D position error due to multipath) with figure 
3 (single channel multipath error). 

The separate channel tracking approach as presented in this paper 
compared to the simultaneous (sub)optimal Kalman filter may have a distinct 
implementation advantage in case the GPS receiver has less than four 
tracking channels. The problem of 'asynchronized filter feeding' can be 
solved by separate channel filtering, fitting a nth order function through 
the filter output data (e.g. a second-order function as described in 
Chapter 5.D.2.2 and finally executing a position flash computation 
(Chapter 5.C) or even using a simultaneous filter using synchronized data 
derived from the fitting function. 
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